module ice_comp_mct

!---------------------------------------------------------------------------
!BOP
!
! !MODULE: ice_comp_mct
!
! !DESCRIPTION:
! CICE interface routine for the ccsm cpl7 mct system
!
! !USES:

  use shr_kind_mod, only : r8 => shr_kind_r8
  use shr_sys_mod,  only : shr_sys_abort, shr_sys_flush
! use shr_mem_mod,  only : shr_get_memusage, shr_init_memusage
  use shr_file_mod, only : shr_file_getlogunit, shr_file_getloglevel,  &
		           shr_file_setloglevel, shr_file_setlogunit
  use mct_mod
#ifdef USE_ESMF_LIB
  use esmf_mod
#else
  use esmf_mod, only: ESMF_clock
#endif

  use seq_flds_mod
  use seq_cdata_mod,   only : seq_cdata, seq_cdata_setptrs
  use seq_infodata_mod,only : seq_infodata_type, seq_infodata_getdata,       &
		              seq_infodata_putdata, seq_infodata_start_type_cont, &
		              seq_infodata_start_type_brnch, seq_infodata_start_type_start
  use seq_timemgr_mod, only : seq_timemgr_eclockgetdata, &
                              seq_timemgr_restartalarmison, &
		              seq_timemgr_eclockdateinsync, &
                              seq_timemgr_stopalarmison
  use perf_mod,        only : t_startf, t_stopf

  use ice_cpl_indices
  use ice_flux,        only : strairxt, strairyt, strocnxt, strocnyt,    &
			      alvdr, alidr, alvdf, alidf, Tref, Qref, Uref, &
                              flat, fsens, flwout, evap, fswabs, fhocn, &
                              fswthru,     &
		              fresh, fsalt, zlvl, uatm, vatm, potT, Tair, Qa,  &
		              rhoa, swvdr, swvdf, swidr, swidf, flw, frain,    &
		              fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt,&
		              sss, tf, wind, fsw, init_flux_atm, init_flux_ocn,&
                              faero
  use ice_state,       only : vice, vsno, aice, trcr, filename_aero, filename_iage, &
                              filename_volpn, filename_FY, filename_lvl, &
                              tr_aero, tr_iage, tr_FY, tr_pond, tr_lvl
  use ice_domain_size, only : nx_global, ny_global, block_size_x, block_size_y, max_blocks
  use ice_domain,      only : nblocks, blocks_ice, halo_info, distrb_info, profile_barrier
  use ice_blocks,      only : block, get_block, nx_block, ny_block
  use ice_grid,        only : tlon, tlat, tarea, tmask, anglet, hm, ocn_gridcell_frac, &
 		              grid_type, t2ugrid_vector
  use ice_constants,   only : c0, c1, puny, tffresh, spval_dbl, rad_to_deg, radius, &
		              field_loc_center, field_type_scalar, field_type_vector, c100
  use ice_communicate, only : my_task, master_task, lprint_stats, MPI_COMM_ICE
  use ice_calendar,    only : idate, mday, time, month, daycal, secday, &
		              sec, dt, dt_dyn, xndt_dyn, calendar,      &
                              calendar_type, nextsw_cday, days_per_year,&
                              get_daycal, leap_year_count
  use ice_orbital,     only : eccen, obliqr, lambm0, mvelpp
  use ice_timers
  use ice_probability, only : init_numIceCells, print_numIceCells,  &
 			      write_numIceCells, accum_numIceCells2

  use ice_kinds_mod,   only : int_kind, dbl_kind, char_len_long, log_kind
!  use ice_init
  use ice_boundary,    only : ice_HaloUpdate 
  use ice_scam,        only : scmlat, scmlon, single_column
  use ice_fileunits,   only : nu_diag
  use ice_dyn_evp,     only :  kdyn
  use ice_prescribed_mod
  use ice_step_mod
  use CICE_RunMod
  use ice_global_reductions
  use ice_broadcast

! !PUBLIC MEMBER FUNCTIONS:
  implicit none
  public :: ice_init_mct
  public :: ice_run_mct
  public :: ice_final_mct
  SAVE
  private                              ! By default make data private
!
! ! PUBLIC DATA:
!
! !REVISION HISTORY:
! Author: Jacob Sewall, Mariana Vertenstein
!
!EOP
! !PRIVATE MEMBER FUNCTIONS:
  private :: ice_export_mct
  private :: ice_import_mct
  private :: ice_SetGSMap_mct
  private :: ice_domain_mct

!
! !PRIVATE VARIABLES

  integer (kind=int_kind) :: ICEID       

  logical (kind=log_kind) :: atm_aero

!=======================================================================

contains

!=======================================================================
!BOP
!
! !IROUTINE: ice_init_mct
!
! !INTERFACE:
  subroutine ice_init_mct( EClock, cdata_i, x2i_i, i2x_i, NLFilename )
!
! !DESCRIPTION:
! Initialize thermodynamic ice model and obtain relevant atmospheric model
! arrays back from driver 
!
! !USES:

    use CICE_InitMod
    use ice_restart, only: runid, runtype, restart_dir, restart_format
    use ice_history, only: history_dir, history_file
!
! !ARGUMENTS:
    type(ESMF_Clock)         , intent(in)    :: EClock
    type(seq_cdata)          , intent(inout) :: cdata_i
    type(mct_aVect)          , intent(inout) :: x2i_i, i2x_i
    character(len=*), optional  , intent(in) :: NLFilename ! Namelist filename
!
! !LOCAL VARIABLES:
!
    type(mct_gsMap)             , pointer :: gsMap_ice
    type(mct_gGrid)             , pointer :: dom_i
    type(seq_infodata_type)     , pointer :: infodata   ! Input init object
    integer                               :: lsize

    character(len=256) :: drvarchdir         ! driver archive directory
    character(len=32)  :: starttype          ! infodata start type
    integer            :: start_ymd          ! Start date (YYYYMMDD)
    integer            :: start_tod          ! start time of day (s)
    integer            :: ref_ymd            ! Reference date (YYYYMMDD)
    integer            :: ref_tod            ! reference time of day (s)
    integer            :: iyear              ! yyyy
    integer            :: dtime              ! time step
    integer            :: shrlogunit,shrloglev ! old values
    integer            :: iam,ierr
    integer            :: lbnum
    integer            :: daycal(13)  !number of cumulative days per month
    integer            :: nleaps      ! number of leap days before current year
    integer            :: mpicom_loc  ! temporary mpicom

    real(r8) :: mrss, mrss0,msize,msize0

! !REVISION HISTORY:
! Author: Mariana Vertenstein
!EOP
!-----------------------------------------------------------------------

    !--------------------------------------------------------------------------
    ! Determine attribute vector indices
    !--------------------------------------------------------------------------

    call ice_cpl_indices_set()

    !---------------------------------------------------------------------------
    ! Set cdata pointers
    !---------------------------------------------------------------------------

    call seq_cdata_setptrs(cdata_i, ID=ICEID, mpicom=mpicom_loc, &
         gsMap=gsMap_ice, dom=dom_i, infodata=infodata)

    ! Determine time of next atmospheric shortwave calculation
    call seq_infodata_GetData(infodata, nextsw_cday=nextsw_cday )

    ! Determine if aerosols are coming from the coupler
    call seq_infodata_GetData(infodata, atm_aero=atm_aero )

    ! Determine orbital parameters
    call seq_infodata_GetData(infodata, orb_eccen=eccen, orb_mvelpp=mvelpp, &
                              orb_lambm0=lambm0, orb_obliqr=obliqr)

    !   call shr_init_memusage()

    !---------------------------------------------------------------------------
    ! use infodata to determine type of run
    !---------------------------------------------------------------------------

    ! Preset single column values

    single_column = .false.
    scmlat = -999.
    scmlon = -999.

    call seq_infodata_GetData( infodata, case_name=runid   ,  &  
       single_column=single_column ,scmlat=scmlat,scmlon=scmlon)
    call seq_infodata_GetData( infodata, start_type=starttype)

    if (     trim(starttype) == trim(seq_infodata_start_type_start)) then
       runtype = "initial"
    else if (trim(starttype) == trim(seq_infodata_start_type_cont) ) then
       runtype = "continue"
    else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then
       runtype = "branch"
    else
       write(nu_diag,*) 'ice_comp_mct ERROR: unknown starttype'
       call shr_sys_abort()
    end if

    ! Set nextsw_cday to -1 for continue and branch runs.

    if (trim(runtype) /= 'initial') nextsw_cday = -1

    !=============================================================
    ! Set ice dtime to ice coupling frequency
    !=============================================================

    call seq_timemgr_EClockGetData(EClock, dtime=dtime, calendar=calendar_type)
    dt = real(dtime)

    !=============================================================
    ! Initialize cice because grid information is needed for
    ! creation of GSMap_ice.  cice_init also sets time manager info
    !=============================================================

    call t_startf ('cice_init')
    call cice_init( mpicom_loc )
    call t_stopf ('cice_init')
    call init_numIceCells

    !---------------------------------------------------------------------------
    ! Reset shr logging to my log file
    !---------------------------------------------------------------------------

    call shr_file_getLogUnit (shrlogunit)
    call shr_file_getLogLevel(shrloglev)
    call shr_file_setLogUnit (nu_diag)
   
    !---------------------------------------------------------------------------
    ! use EClock to reset calendar information on initial start
    !---------------------------------------------------------------------------

    ! - the following logic duplicates the logic for the concurrent system - 
    ! cice_init is called then init_cpl is called where the start date is received
    ! from the flux coupler
    ! - in the following calculation for the variable time, iyear is used
    ! rather than (iyear-1) (as in init_cpl) since the sequential system permits you 
    ! to start with year "0" not year "1"
    ! - on restart run 
    !   - istep0, time and time_forc are read from restart file
    !   - istep1 is set to istep0
    !   - idate is determined from time via the call to calendar (see below)
    ! - on initial run 
    !   - iyear, month and mday obtained from sync clock
    !   - time determined from iyear, month and mday
    !   - istep0 and istep1 are set to 0 

    if (runtype == 'initial') then
       call seq_timemgr_EClockGetData(EClock, &
            start_ymd=start_ymd, start_tod=start_tod,       &
            ref_ymd=ref_ymd, ref_tod=ref_tod)

       if (ref_ymd /= start_ymd .or. ref_tod /= start_tod) then
          if (my_task == master_task) then
             write(nu_diag,*) 'ice_comp_mct: ref_ymd ',ref_ymd, &
                  ' must equal start_ymd ',start_ymd
             write(nu_diag,*) 'ice_comp_mct: ref_ymd ',ref_tod, &
                  ' must equal start_ymd ',start_tod
          end if
          call shr_sys_abort()
       end if

       if (my_task == master_task) then
          write(nu_diag,*) '(ice_init_mct) idate from sync clock = ', &
               start_ymd
          write(nu_diag,*) '(ice_init_mct)   tod from sync clock = ', &
               start_tod
          write(nu_diag,*) &
               '(ice_init_mct) resetting idate to match sync clock'
       end if

       idate = start_ymd
       iyear = (idate/10000)                     ! integer year of basedate
       month = (idate-iyear*10000)/100           ! integer month of basedate
       mday  =  idate-iyear*10000-month*100-1    ! day of month of basedate
                                                 ! (starts at 0)

       call get_daycal(year=iyear,days_per_year_in=days_per_year, &
                       daycal_out=daycal)

       nleaps = leap_year_count(iyear)   ! this sets nleaps in ice_calendar
       time  = (((iyear)*days_per_year  + nleaps + daycal(month)+mday)*secday) &
             + start_tod    

       call shr_sys_flush(nu_diag)
    end if

    call calendar(time)     ! update calendar info
 
    !---------------------------------------------------------------------------
    ! Initialize MCT attribute vectors and indices
    !---------------------------------------------------------------------------

    call t_startf ('cice_mct_init')

    ! Initialize ice gsMap

    call ice_SetGSMap_mct( MPI_COMM_ICE, ICEID, GSMap_ice ) 	
    lsize = mct_gsMap_lsize(gsMap_ice, MPI_COMM_ICE)

    ! Initialize mct ice domain (needs ice initialization info)

    call ice_domain_mct( lsize, gsMap_ice, dom_i )

    ! Inialize mct attribute vectors

    call mct_aVect_init(x2i_i, rList=seq_flds_x2i_fields, lsize=lsize)
    call mct_aVect_zero(x2i_i)

    call mct_aVect_init(i2x_i, rList=seq_flds_i2x_fields, lsize=lsize) 
    call mct_aVect_zero(i2x_i)

    !-----------------------------------------------------------------
    ! Prescribed ice initialization
    !-----------------------------------------------------------------

    call ice_prescribed_init(ICEID, gsmap_ice, dom_i)

    !-----------------------------------------------------------------
    ! Get ready for coupling
    !-----------------------------------------------------------------

    call coupling_prep

    !---------------------------------------------------------------------------
    ! Fill in export state for driver
    !---------------------------------------------------------------------------

    call ice_export_mct (i2x_i)  !Send initial state to driver
    call seq_infodata_PutData( infodata, ice_prognostic=.true., &
      ice_nx = nx_global, ice_ny = ny_global )
    call t_stopf ('cice_mct_init')

    !---------------------------------------------------------------------------
    ! Reset shr logging to original values
    !---------------------------------------------------------------------------

    call shr_file_setLogUnit (shrlogunit)
    call shr_file_setLogLevel(shrloglev)

    call ice_timer_stop(timer_total) ! time entire run
!   call shr_get_memusage(msize,mrss)
!   call shr_mpi_max(mrss, mrss0, MPI_COMM_ICE,'ice_init_mct mrss0')
!   call shr_mpi_max(msize,msize0,MPI_COMM_ICE,'ice_init_mct msize0')
!   if(my_task == 0) then
!     write(shrlogunit,105) 'ice_init_mct: memory_write: model date = ',start_ymd,start_tod, &
!           ' memory = ',msize0,' MB (highwater)    ',mrss0,' MB (usage)'
!   endif
 
  105  format( A, 2i8, A, f10.2, A, f10.2, A)

  end subroutine ice_init_mct

!---------------------------------------------------------------------------
!BOP
!
! !IROUTINE: ice_run_mct
!
! !INTERFACE:
  subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i )
!
! !DESCRIPTION:
! Run thermodynamic CICE
!
! !USES:
    use ice_history
    use ice_restart
    use ice_diagnostics
    use ice_aerosol, only: write_restart_aero
    use ice_age, only: write_restart_age
    use ice_meltpond, only: write_restart_pond
    use ice_FY, only: write_restart_FY
    use ice_lvl, only: write_restart_lvl
    use ice_restoring, only: restore_ice, ice_HaloRestore
    use ice_shortwave, only: init_shortwave

! !ARGUMENTS:
    type(ESMF_Clock),intent(in)    :: EClock
    type(seq_cdata), intent(inout) :: cdata_i
    type(mct_aVect), intent(inout) :: x2i_i
    type(mct_aVect), intent(inout) :: i2x_i

! !LOCAL VARIABLES:
    integer :: k             ! index
    logical :: rstwr         ! .true. ==> write a restart file
    logical :: stop_now      ! .true. ==> stop at the end of this run phase
    integer :: ymd           ! Current date (YYYYMMDD)
    integer :: tod           ! Current time of day (sec)
    integer :: yr_sync       ! Sync current year
    integer :: mon_sync      ! Sync current month
    integer :: day_sync      ! Sync current day
    integer :: tod_sync      ! Sync current time of day (sec)
    integer :: ymd_sync      ! Current year of sync clock
    integer :: shrlogunit,shrloglev ! old values
    integer :: lbnum
    integer :: n
    type(mct_gGrid)        , pointer :: dom_i
    type(seq_infodata_type), pointer :: infodata   
    type(mct_gsMap)        , pointer :: gsMap_i
    character(len=char_len_long) :: fname
    character(len=char_len_long) :: string1, string2
    character(len=*), parameter  :: SubName = "ice_run_mct"

    real(r8) :: mrss, mrss0,msize,msize0
    logical, save :: first_time = .true.

!
! !REVISION HISTORY:
! Author: Jacob Sewall, Mariana Vertenstein
!
!EOP
!---------------------------------------------------------------------------

    call ice_timer_start(timer_total) ! time entire run

    !---------------------------------------------------------------------------
    ! Reset shr logging to my log file
    !---------------------------------------------------------------------------

    call shr_file_getLogUnit (shrlogunit)
    call shr_file_getLogLevel(shrloglev)
    call shr_file_setLogUnit (nu_diag)
   
    call seq_cdata_setptrs(cdata_i, infodata=infodata, dom=dom_i, &
         gsMap=gsMap_i)

    ! Determine time of next atmospheric shortwave calculation
    call seq_infodata_GetData(infodata, nextsw_cday=nextsw_cday )

    ! Determine orbital parameters
    call seq_infodata_GetData(infodata, orb_eccen=eccen, orb_mvelpp=mvelpp, &
                              orb_lambm0=lambm0, orb_obliqr=obliqr)

    !-------------------------------------------------------------------
    ! get import state
    !-------------------------------------------------------------------
    
    call t_startf ('cice_import')
    call ice_timer_start(timer_cplrecv)
    call ice_import_mct( x2i_i )
    call ice_timer_stop(timer_cplrecv)
    call t_stopf ('cice_import')
 
    !--------------------------------------------------------------------
    ! timestep update
    !--------------------------------------------------------------------

    call ice_timer_start(timer_step)

    istep  = istep  + 1    ! update time step counters
    istep1 = istep1 + 1
    time = time + dt       ! determine the time and date
    call calendar(time)    ! at the end of the timestep

    if (resttype == 'old' .and. istep == 1) &
       call init_shortwave    ! initialize radiative transfer using current swdn

    !-----------------------------------------------------------------
    ! restoring on grid boundaries
    !-----------------------------------------------------------------

    if (restore_ice) call ice_HaloRestore

    call t_startf ('cice_initmd')
    call init_mass_diags   ! diagnostics per timestep
    call t_stopf ('cice_initmd')

    if(prescribed_ice) then  ! read prescribed ice
       call t_startf ('cice_presc')
       call ice_prescribed_run(idate, sec)
       call t_stopf ('cice_presc')
    endif
    
    call init_flux_atm        ! initialize atmosphere fluxes sent to coupler
    call init_flux_ocn        ! initialize ocean fluxes sent to coupler

    !-----------------------------------------------------------------
    ! Scale radiation fields
    !-----------------------------------------------------------------

    call t_startf ('cice_prep_radiation')
    call ice_timer_start(timer_sw)
    call prep_radiation(dt)
    call ice_timer_stop(timer_sw)
    call t_stopf ('cice_prep_radiation')
    
    !-----------------------------------------------------------------
    ! thermodynamics1
    !-----------------------------------------------------------------

    call t_startf ('cice_therm1')
    call step_therm1(dt)
    call t_stopf ('cice_therm1')
    
    !-----------------------------------------------------------------
    ! thermodynamics2
    !-----------------------------------------------------------------

    if (.not.prescribed_ice) then
       call t_startf ('cice_therm2')
       call step_therm2 (dt)  ! post-coupler thermodynamics
       call t_stopf ('cice_therm2')
    end if

   !-----------------------------------------------------------------
   ! dynamics, transport, ridging
   !-----------------------------------------------------------------

    if (.not.prescribed_ice .and. kdyn>0) then
       if (xndt_dyn > c1) then
          call t_startf ('cice_dyn')
          do k = 1, nint(xndt_dyn)
             call step_dynamics(dt_dyn,dt) ! dynamics, transport, ridging
          enddo
          call t_stopf ('cice_dyn')
       else
          if (mod(time, dt_dyn) == c0) then
             call t_startf ('cice_dyn')
             call step_dynamics(dt_dyn,dt) ! dynamics, transport, ridging
             call t_stopf ('cice_dyn')
          endif
       endif
    endif ! not prescribed_ice
    
    !-----------------------------------------------------------------
    ! radiation
    !-----------------------------------------------------------------

    call t_startf ('cice_radiation')
    call ice_timer_start(timer_sw)
    call step_radiation(dt)
    call ice_timer_stop(timer_sw)
    call t_stopf ('cice_radiation')
    
    !-----------------------------------------------------------------
    ! get ready for coupling
    !-----------------------------------------------------------------

    call coupling_prep

    call ice_timer_stop(timer_step)

    !-----------------------------------------------------------------
    ! write data
    !-----------------------------------------------------------------
    
    call t_startf ('cice_diag')
    call ice_timer_start(timer_diags)
    if (mod(istep,diagfreq) == 0) call runtime_diags(dt) ! log file
    call ice_timer_stop(timer_diags)
    call t_stopf ('cice_diag')
    
    call t_startf ('cice_hist')
    call ice_timer_start(timer_hist)
#if (defined _NOIO)
!  Not enought memory on BGL to write a history file yet! 
!    call ice_write_hist (dt)    ! history file
#else
    call ice_write_hist (dt)    ! history file
#endif
    call ice_timer_stop(timer_hist)
    call t_stopf ('cice_hist')
 
    !--------------------------------------------
    ! Accumualate the number of active ice cells
    !--------------------------------------------
    call accum_numIceCells2(aice)

    rstwr = seq_timemgr_RestartAlarmIsOn(EClock)
    if (rstwr) then
       call seq_timemgr_EClockGetData(EClock, curr_ymd=ymd_sync, curr_tod=tod_sync, &
          curr_yr=yr_sync,curr_mon=mon_sync,curr_day=day_sync)
       fname = restart_filename(yr_sync, mon_sync, day_sync, tod_sync)

       if (my_task == master_task) then
          write(nu_diag,*) &
            'ice_comp_mct: calling dumpfile for restart filename= ', fname
       endif

       if (restart_format /= 'nc') then

          if (tr_pond) then
               n = index(fname,'cice.r') + 6
               string1 = trim(fname(1:n-1))
               string2 = trim(fname(n:lenstr(fname)))
               write(filename_volpn,'(a,a,a,a)') &
                  string1(1:lenstr(string1)),'.volpn', &
                  string2(1:lenstr(string2))
          endif
          if (tr_aero) then
               n = index(fname,'cice.r') + 6
               string1 = trim(fname(1:n-1))
               string2 = trim(fname(n:lenstr(fname)))
               write(filename_aero,'(a,a,a,a)') &
                  string1(1:lenstr(string1)),'.aero', &
                  string2(1:lenstr(string2))
          endif
          if (tr_iage) then
               n = index(fname,'cice.r') + 6
               string1 = trim(fname(1:n-1))
               string2 = trim(fname(n:lenstr(fname)))
               write(filename_iage,'(a,a,a,a)') &
                  string1(1:lenstr(string1)),'.age', &
                  string2(1:lenstr(string2))
          endif
          if (tr_FY) then
               n = index(fname,'cice.r') + 6
               string1 = trim(fname(1:n-1))
               string2 = trim(fname(n:lenstr(fname)))
               write(filename_FY,'(a,a,a,a)') &
                  string1(1:lenstr(string1)),'.FY', &
                  string2(1:lenstr(string2))
          endif
          if (tr_lvl) then
               n = index(fname,'cice.r') + 6
               string1 = trim(fname(1:n-1))
               string2 = trim(fname(n:lenstr(fname)))
               write(filename_lvl,'(a,a,a,a)') &
                  string1(1:lenstr(string1)),'.lvl', &
                  string2(1:lenstr(string2))
          endif

       endif ! restart_format

#if (defined _NOIO)
!  Not enought memory on BGL to call dumpfile  file yet! 
!       call dumpfile(fname)
#else
       call ice_timer_start(timer_readwrite)
       call dumpfile(fname)
       if (restart_format /= 'nc') then
          if (tr_aero) call write_restart_aero(filename_aero)
          if (tr_iage) call write_restart_age(filename_iage)
          if (tr_FY)   call write_restart_FY(filename_FY)
          if (tr_lvl)  call write_restart_lvl(filename_lvl)
          if (tr_pond) call write_restart_pond(filename_volpn)
       endif
       call ice_timer_stop(timer_readwrite)
#endif
    end if

    !-----------------------------------------------------------------
    ! send export state to driver 
    !-----------------------------------------------------------------
    
    call t_startf ('cice_export')
    call ice_timer_start(timer_cplsend)
    call ice_export_mct ( i2x_i )
    call ice_timer_stop(timer_cplsend)
    call t_stopf ('cice_export')
    
    !--------------------------------------------------------------------
    ! check that internal clock is in sync with master clock
    !--------------------------------------------------------------------

    tod = sec
    ymd = idate
    if (.not. seq_timemgr_EClockDateInSync( EClock, ymd, tod )) then
       call seq_timemgr_EClockGetData( EClock, curr_ymd=ymd_sync, &
          curr_tod=tod_sync )
       write(nu_diag,*)' cice ymd=',ymd     ,'  cice tod= ',tod
       write(nu_diag,*)' sync ymd=',ymd_sync,'  sync tod= ',tod_sync
       call shr_sys_abort( SubName// &
          ":: Internal sea-ice clock not in sync with Sync Clock")
    end if
   
    ! reset shr logging to my original values

    call shr_file_setLogUnit (shrlogunit)
    call shr_file_setLogLevel(shrloglev)

    !-------------------------------------------------------------------
    ! stop timers and print timer info
    !-------------------------------------------------------------------
    ! Need to have this logic here instead of in ice_final_mct since 
    ! the ice_final_mct.F90 will still be called even in aqua-planet mode
    ! Could put this logic in the driver - but it seems easier here 

    ! Need to stop this at the end of every run phase in a coupled run.
    call ice_timer_stop(timer_total)        ! stop timing

    stop_now = seq_timemgr_StopAlarmIsOn( EClock )
    if (stop_now) then
       call ice_timer_print_all(stats=.true.) ! print timing information
       if(lprint_stats) then 
          call write_numIceCells()
       endif
       call release_all_fileunits
    end if
    
!   if(tod == 0) then
!      call shr_get_memusage(msize,mrss)
!      call shr_mpi_max(mrss, mrss0, MPI_COMM_ICE,'ice_run_mct mrss0')
!      call shr_mpi_max(msize,msize0,MPI_COMM_ICE,'ice_run_mct msize0')
!      if(my_task == 0 ) then
!          write(shrlogunit,105) 'ice_run_mct: memory_write: model date = ',ymd,tod, &
!               ' memory = ',msize0,' MB (highwater)    ',mrss0,' MB (usage)'
!      endif
!   endif
 
  105  format( A, 2i8, A, f10.2, A, f10.2, A)

  end subroutine ice_run_mct

!---------------------------------------------------------------------------
!BOP
!
! !IROUTINE: ice_final_mct
!
! !INTERFACE:
  subroutine ice_final_mct( )
!
! !DESCRIPTION:
! Finalize CICE
!
! !USES:
!
!------------------------------------------------------------------------------
!BOP
!
! !ARGUMENTS:
!
! !REVISION HISTORY:
!
!EOP
!---------------------------------------------------------------------------

  end subroutine ice_final_mct

!===============================================================================

  subroutine ice_SetGSMap_mct( mpicom, ID, gsMap_ice )

    !-------------------------------------------------------------------
    !
    ! Arguments
    !
    implicit none
    integer        , intent(in)    :: mpicom
    integer        , intent(in)    :: ID
    type(mct_gsMap), intent(inout) :: gsMap_ice
    !
    ! Local variables
    !
    integer,allocatable :: gindex(:)
    integer     :: lat
    integer     :: lon
    integer     :: i, j, iblk, n, gi
    integer     :: lsize,gsize
    integer     :: ier
    integer     :: ilo, ihi, jlo, jhi ! beginning and end of physical domain
    type(block) :: this_block         ! block information for current block
    !-------------------------------------------------------------------

    ! Build the CICE grid numbering for MCT
    ! NOTE:  Numbering scheme is: West to East and South to North
    ! starting at south pole.  Should be the same as what's used
    ! in SCRIP

    ! number the local grid

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
          do i = ilo, ihi
             n = n+1
          enddo !i
       enddo    !j
    enddo        !iblk
    lsize = n

! not valid for padded decomps
!    lsize = block_size_x*block_size_y*nblocks
    gsize = nx_global*ny_global

    allocate(gindex(lsize),stat=ier)
    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
          do i = ilo, ihi
             n = n+1
             lon = this_block%i_glob(i)
             lat = this_block%j_glob(j)
             gi = (lat-1)*nx_global + lon
             gindex(n) = gi
          enddo !i
       enddo    !j
    enddo        !iblk
    
    call mct_gsMap_init( gsMap_ice, gindex, mpicom, ID, lsize, gsize )

    deallocate(gindex)

  end subroutine ice_SetGSMap_mct


!===============================================================================

  subroutine ice_export_mct( i2x_i )   

    !-----------------------------------------------------
    type(mct_aVect)   , intent(inout) :: i2x_i

    integer :: i, j, iblk, n, ij 
    integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain
    integer (kind=int_kind)                                :: icells ! number of ocean/ice cells
    integer (kind=int_kind), dimension (nx_block*ny_block) :: indxi  ! compressed indices in i
    integer (kind=int_kind), dimension (nx_block*ny_block) :: indxj  ! compressed indices in i

    real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
        Tsrf  &      ! surface temperature
     ,  tauxa &      ! atmo/ice stress
     ,  tauya &
     ,  tauxo &      ! ice/ocean stress
     ,  tauyo &
     ,  ailohi       ! fractional ice area

    real (kind=dbl_kind) :: &
       workx, worky           ! tmps for converting grid

    type(block)        :: this_block                           ! block information for current block
    logical :: flag
    !-----------------------------------------------------

    flag=.false.

    !calculate ice thickness from aice and vice. Also
    !create Tsrf from the first tracer (trcr) in ice_state.F

    !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky)
    do iblk = 1, nblocks
       do j = 1, ny_block
       do i = 1, nx_block
             
          ! ice fraction
          ailohi(i,j,iblk) = min(aice(i,j,iblk), c1)

          ! surface temperature
          Tsrf(i,j,iblk)  = Tffresh + trcr(i,j,1,iblk)             !Kelvin (original ???)
          
          ! wind stress  (on POP T-grid:  convert to lat-lon)
          workx = strairxT(i,j,iblk)                             ! N/m^2
          worky = strairyT(i,j,iblk)                             ! N/m^2
          tauxa(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) &
                          - worky*sin(ANGLET(i,j,iblk))
          tauya(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
                          + workx*sin(ANGLET(i,j,iblk))
          
          ! ice/ocean stress (on POP T-grid:  convert to lat-lon)
          workx = -strocnxT(i,j,iblk)                            ! N/m^2
          worky = -strocnyT(i,j,iblk)                            ! N/m^2
          tauxo(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) &
                          - worky*sin(ANGLET(i,j,iblk))
          tauyo(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
                          + workx*sin(ANGLET(i,j,iblk))

       enddo
       enddo
    enddo
    !$OMP END PARALLEL DO

    do iblk = 1, nblocks
       do j = 1, ny_block
       do i = 1, nx_block
          if (tmask(i,j,iblk) .and. ailohi(i,j,iblk) < c0 ) then
             flag = .true.
          endif
       end do
       end do
    end do
    if (flag) then
       do iblk = 1, nblocks
          do j = 1, ny_block
          do i = 1, nx_block
             if (tmask(i,j,iblk) .and. ailohi(i,j,iblk) < c0 ) then
                write(nu_diag,*) &
                     ' (ice) send: ERROR ailohi < 0.0 ',i,j,ailohi(i,j,iblk)
                call shr_sys_flush(nu_diag)
             endif
          end do
          end do
       end do
    endif

    ! Fill export state i2x_i

     i2x_i%rAttr(:,:) = spval_dbl

     n=0
     do iblk = 1, nblocks
         this_block = get_block(blocks_ice(iblk),iblk)         
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
         jhi = this_block%jhi

         do j = jlo, jhi
         do i = ilo, ihi

            n = n+1

            !-------states-------------------- 
            i2x_i%rAttr(index_i2x_Si_ifrac ,n)    = ailohi(i,j,iblk)   

            if ( tmask(i,j,iblk) .and. ailohi(i,j,iblk) > c0 ) then
               !-------states-------------------- 
               i2x_i%rAttr(index_i2x_Si_t     ,n)    = Tsrf(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_avsdr ,n)    = alvdr(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_anidr ,n)    = alidr(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_avsdf ,n)    = alvdf(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_anidf ,n)    = alidf(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_u10  ,n)     = Uref(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_tref  ,n)    = Tref(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_qref  ,n)    = Qref(i,j,iblk)
               i2x_i%rAttr(index_i2x_Si_snowh ,n)    = vsno(i,j,iblk) &
                                                     / ailohi(i,j,iblk)
            
               !--- a/i fluxes computed by ice
               i2x_i%rAttr(index_i2x_Faii_taux ,n)   = tauxa(i,j,iblk)    
               i2x_i%rAttr(index_i2x_Faii_tauy ,n)   = tauya(i,j,iblk)    
               i2x_i%rAttr(index_i2x_Faii_lat  ,n)   = flat(i,j,iblk)     
               i2x_i%rAttr(index_i2x_Faii_sen  ,n)   = fsens(i,j,iblk)    
               i2x_i%rAttr(index_i2x_Faii_lwup ,n)   = flwout(i,j,iblk)   
               i2x_i%rAttr(index_i2x_Faii_evap ,n)   = evap(i,j,iblk)     
               i2x_i%rAttr(index_i2x_Faii_swnet,n)   = fswabs(i,j,iblk)
            
               !--- i/o fluxes computed by ice
               i2x_i%rAttr(index_i2x_Fioi_melth,n)   = fhocn(i,j,iblk)
               i2x_i%rAttr(index_i2x_Fioi_swpen,n)   = fswthru(i,j,iblk) ! hf from melting          
               i2x_i%rAttr(index_i2x_Fioi_meltw,n)   = fresh(i,j,iblk)   ! h2o flux from melting    ???
               i2x_i%rAttr(index_i2x_Fioi_salt ,n)   = fsalt(i,j,iblk)   ! salt flux from melting   ???
               i2x_i%rAttr(index_i2x_Fioi_taux ,n)   = tauxo(i,j,iblk)   ! stress : i/o zonal       ???
               i2x_i%rAttr(index_i2x_Fioi_tauy ,n)   = tauyo(i,j,iblk)   ! stress : i/o meridional  ???
            end if
         enddo    !i
         enddo    !j
     enddo        !iblk

  end subroutine ice_export_mct

!==============================================================================

  subroutine ice_import_mct( x2i_i )

    !-----------------------------------------------------

    implicit none
    type(mct_aVect)   , intent(inout) :: x2i_i

    integer :: i, j, iblk, n
    integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain
    type(block) :: this_block      ! block information for current block

    integer,parameter :: nflds=15,nfldv=6
    real (kind=dbl_kind),allocatable :: aflds(:,:,:,:)

    real (kind=dbl_kind) :: &
         workx, worky
    logical (kind=log_kind) :: &
         first_call = .true.
    logical (kind=log_kind) :: &
         prescribed_aero_tmp

    !-----------------------------------------------------

    ! Note that the precipitation fluxes received  from the coupler
    ! are in units of kg/s/m^2 which is what CICE requires.
    ! Note also that the read in below includes only values needed
    ! by the thermodynamic component of CICE.  Variables uocn, vocn,
    ! ss_tltx, and ss_tlty are excluded. Also, because the SOM and
    ! DOM don't  compute SSS.   SSS is not read in and is left at
    ! the initilized value (see ice_flux.F init_coupler_flux) of
    ! 34 ppt

    ! Use aflds to gather the halo updates of multiple fields
    ! Need to separate the scalar from the vector halo updates

    allocate(aflds(nx_block,ny_block,nflds,nblocks))
    aflds = c0

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi

       do j = jlo, jhi
       do i = ilo, ihi

          n = n+1
          aflds(i,j, 1,iblk)   = x2i_i%rAttr(index_x2i_So_t,n)
          aflds(i,j, 2,iblk)   = x2i_i%rAttr(index_x2i_So_s,n)
          aflds(i,j, 3,iblk)   = x2i_i%rAttr(index_x2i_Sa_z,n)
          aflds(i,j, 4,iblk)   = x2i_i%rAttr(index_x2i_Sa_ptem,n)
          aflds(i,j, 5,iblk)   = x2i_i%rAttr(index_x2i_Sa_tbot,n)
          aflds(i,j, 6,iblk)   = x2i_i%rAttr(index_x2i_Sa_shum,n)
          aflds(i,j, 7,iblk)   = x2i_i%rAttr(index_x2i_Sa_dens,n)
          aflds(i,j, 8,iblk)   = x2i_i%rAttr(index_x2i_Fioo_q,n)
          aflds(i,j, 9,iblk)   = x2i_i%rAttr(index_x2i_Faxa_swvdr,n)
          aflds(i,j,10,iblk)   = x2i_i%rAttr(index_x2i_Faxa_swndr,n)
          aflds(i,j,11,iblk)   = x2i_i%rAttr(index_x2i_Faxa_swvdf,n)
          aflds(i,j,12,iblk)   = x2i_i%rAttr(index_x2i_Faxa_swndf,n)
          aflds(i,j,13,iblk)   = x2i_i%rAttr(index_x2i_Faxa_lwdn,n)
          aflds(i,j,14,iblk)   = x2i_i%rAttr(index_x2i_Faxa_rain,n)
          aflds(i,j,15,iblk)   = x2i_i%rAttr(index_x2i_Faxa_snow,n)

         enddo    !i
         enddo    !j

     enddo        !iblk

     if (.not.prescribed_ice) then
        call t_startf ('cice_imp_halo')
        call ice_HaloUpdate(aflds, halo_info, field_loc_center, &
                                              field_type_scalar)
        call t_stopf ('cice_imp_halo')
     endif
 
     !$OMP PARALLEL DO PRIVATE(iblk,i,j)
     do iblk = 1, nblocks
        do j = 1,ny_block
        do i = 1,nx_block
           sst  (i,j,iblk)   = aflds(i,j, 1,iblk)
           sss  (i,j,iblk)   = aflds(i,j, 2,iblk)
           zlvl (i,j,iblk)   = aflds(i,j, 3,iblk)
           potT (i,j,iblk)   = aflds(i,j, 4,iblk)
           Tair (i,j,iblk)   = aflds(i,j, 5,iblk)
           Qa   (i,j,iblk)   = aflds(i,j, 6,iblk)
           rhoa (i,j,iblk)   = aflds(i,j, 7,iblk)
           frzmlt (i,j,iblk) = aflds(i,j, 8,iblk)
           swvdr(i,j,iblk)   = aflds(i,j, 9,iblk)
           swidr(i,j,iblk)   = aflds(i,j,10,iblk)
           swvdf(i,j,iblk)   = aflds(i,j,11,iblk)
           swidf(i,j,iblk)   = aflds(i,j,12,iblk)
           flw  (i,j,iblk)   = aflds(i,j,13,iblk)
           frain(i,j,iblk)   = aflds(i,j,14,iblk)
           fsnow(i,j,iblk)   = aflds(i,j,15,iblk)
        enddo    !i
        enddo    !j
     enddo        !iblk
     !$OMP END PARALLEL DO
 
     deallocate(aflds)
     allocate(aflds(nx_block,ny_block,nfldv,nblocks))
     aflds = c0

     n=0
     do iblk = 1, nblocks
        this_block = get_block(blocks_ice(iblk),iblk)
        ilo = this_block%ilo
        ihi = this_block%ihi
        jlo = this_block%jlo
        jhi = this_block%jhi
 
        do j = jlo, jhi
        do i = ilo, ihi
           n = n+1
           aflds(i,j, 1,iblk)   = x2i_i%rAttr(index_x2i_So_u,n)
           aflds(i,j, 2,iblk)   = x2i_i%rAttr(index_x2i_So_v,n)
           aflds(i,j, 3,iblk)   = x2i_i%rAttr(index_x2i_Sa_u,n)
           aflds(i,j, 4,iblk)   = x2i_i%rAttr(index_x2i_Sa_v,n)
           aflds(i,j, 5,iblk)   = x2i_i%rAttr(index_x2i_So_dhdx,n)
           aflds(i,j, 6,iblk)   = x2i_i%rAttr(index_x2i_So_dhdy,n)
        enddo
        enddo
     enddo

     if (.not.prescribed_ice) then
        call t_startf ('cice_imp_halo')
        call ice_HaloUpdate(aflds, halo_info, field_loc_center, &
                                             field_type_vector)
        call t_stopf ('cice_imp_halo')
     endif

     !$OMP PARALLEL DO PRIVATE(iblk,i,j)
     do iblk = 1, nblocks
        do j = 1,ny_block
        do i = 1,nx_block
           uocn (i,j,iblk)   = aflds(i,j, 1,iblk)
           vocn (i,j,iblk)   = aflds(i,j, 2,iblk)
           uatm (i,j,iblk)   = aflds(i,j, 3,iblk)
           vatm (i,j,iblk)   = aflds(i,j, 4,iblk)
           ss_tltx(i,j,iblk) = aflds(i,j, 5,iblk)
           ss_tlty(i,j,iblk) = aflds(i,j, 6,iblk)
        enddo    !i
        enddo    !j
     enddo        !iblk
     !$OMP END PARALLEL DO

     deallocate(aflds)

     !-------------------------------------------------------
     ! Set aerosols from coupler 
     !-------------------------------------------------------
     
      if (first_call) then
         if (tr_aero .and. .not. atm_aero) then
            write(nu_diag,*) 'ice_comp_mct ERROR: atm_aero must be set for tr_aero' 
            call shr_sys_abort()
         end if
	 first_call = .false.
      end if

      n=0
      do iblk = 1, nblocks
         this_block = get_block(blocks_ice(iblk),iblk)         
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
         jhi = this_block%jhi

         do j = jlo, jhi
         do i = ilo, ihi

            n = n+1
            faero(i,j,1,iblk) = x2i_i%rAttr(index_x2i_Faxa_bcphodry,n)

            faero(i,j,2,iblk) = x2i_i%rAttr(index_x2i_Faxa_bcphidry,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_bcphiwet,n)
         ! Combine all of the dust into one category
            faero(i,j,3,iblk) = x2i_i%rAttr(index_x2i_Faxa_dstwet1,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstdry1,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstwet2,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstdry2,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstwet3,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstdry3,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstwet4,n) &
                              + x2i_i%rAttr(index_x2i_Faxa_dstdry4,n)

         enddo    !i
         enddo    !j

      enddo        !iblk


     !-----------------------------------------------------------------
     ! rotate zonal/meridional vectors to local coordinates
     ! compute data derived quantities
     !-----------------------------------------------------------------
     
     ! Vector fields come in on T grid, but are oriented geographically
     ! need to rotate to pop-grid FIRST using ANGLET
     ! then interpolate to the U-cell centers  (otherwise we
     ! interpolate across the pole)
     ! use ANGLET which is on the T grid !
     
     call t_startf ('cice_imp_ocn')
     !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky)
     do iblk = 1, nblocks

       do j = 1,ny_block
       do i = 1,nx_block

          ! ocean
          workx      = uocn  (i,j,iblk) ! currents, m/s 
          worky      = vocn  (i,j,iblk)
          uocn(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) & ! convert to POP grid 
                         + worky*sin(ANGLET(i,j,iblk))
          vocn(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
                         - workx*sin(ANGLET(i,j,iblk))

          workx      = ss_tltx  (i,j,iblk)           ! sea sfc tilt, m/m
          worky      = ss_tlty  (i,j,iblk)
          ss_tltx(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) & ! convert to POP grid 
                            + worky*sin(ANGLET(i,j,iblk))
          ss_tlty(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
                            - workx*sin(ANGLET(i,j,iblk))

          sst(i,j,iblk) = sst(i,j,iblk) - Tffresh       ! sea sfc temp (C)
          Tf (i,j,iblk) = -1.8_dbl_kind                 ! hardwired for NCOM
!         Tf (i,j,iblk) = -depressT*sss(i,j,iblk)       ! freezing temp (C)
!         Tf (i,j,iblk) = -depressT*max(sss(i,j,iblk),ice_ref_salinity)

       enddo
       enddo
      enddo
      !$OMP END PARALLEL DO
      call t_stopf ('cice_imp_ocn')

      ! Interpolate ocean dynamics variables from T-cell centers to 
      ! U-cell centers.

      if (.not.prescribed_ice) then
         call t_startf ('cice_imp_t2u')
         call t2ugrid_vector(uocn)
         call t2ugrid_vector(vocn)
         call t2ugrid_vector(ss_tltx)
         call t2ugrid_vector(ss_tlty)
         call t_stopf ('cice_imp_t2u')
      end if

      ! Atmosphere variables are needed in T cell centers in
      ! subroutine stability and are interpolated to the U grid
      ! later as necessary.

      call t_startf ('cice_imp_atm')
      !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky)
      do iblk = 1, nblocks
         do j = 1, ny_block
         do i = 1, nx_block
  
         ! atmosphere
         workx      = uatm(i,j,iblk) ! wind velocity, m/s
         worky      = vatm(i,j,iblk) 
         uatm (i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) & ! convert to POP grid
                         + worky*sin(ANGLET(i,j,iblk))   ! note uatm, vatm, wind
         vatm (i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) & ! are on the T-grid here
                         - workx*sin(ANGLET(i,j,iblk))

         wind (i,j,iblk) = sqrt(uatm(i,j,iblk)**2 + vatm(i,j,iblk)**2)
         fsw  (i,j,iblk) = swvdr(i,j,iblk) + swvdf(i,j,iblk) &
                         + swidr(i,j,iblk) + swidf(i,j,iblk)
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO
      call t_stopf ('cice_imp_atm')

   end subroutine ice_import_mct

!=======================================================================

  subroutine ice_domain_mct( lsize, gsMap_i, dom_i )

    !-------------------------------------------------------------------
    !
    ! Arguments
    !
    integer        , intent(in)    :: lsize
    type(mct_gsMap), intent(in)    :: gsMap_i
    type(mct_ggrid), intent(inout) :: dom_i     
    !
    ! Local Variables
    !
    integer :: i, j, iblk, n, gi           ! indices
    integer :: ilo, ihi, jlo, jhi          ! beginning and end of physical domain
    real(dbl_kind), pointer :: work_dom(:) ! temporary
    real(dbl_kind), pointer :: data(:)     ! temporary
    integer       , pointer :: idata(:)    ! temporary
    type(block)             :: this_block  ! block information for current block
    !-------------------------------------------------------------------
    !
    ! Initialize mct domain type
    ! lat/lon in degrees,  area in radians^2, mask is 1 (ocean), 0 (non-ocean)
    !
    call mct_gGrid_init( GGrid=dom_i, CoordChars=trim(seq_flds_dom_coord), &
       OtherChars=trim(seq_flds_dom_other), lsize=lsize )
    call mct_aVect_zero(dom_i%data)
    !  
    allocate(data(lsize))
    !
    ! Determine global gridpoint number attribute, GlobGridNum, which is set automatically by MCT
    !
    call mct_gsMap_orderedPoints(gsMap_i, my_task, idata)
    call mct_gGrid_importIAttr(dom_i,'GlobGridNum',idata,lsize)
    !
    ! Determine domain (numbering scheme is: West to East and South to North to South pole)
    ! Initialize attribute vector with special value
    !
    data(:) = -9999.0_R8 
    call mct_gGrid_importRAttr(dom_i,"lat"  ,data,lsize) 
    call mct_gGrid_importRAttr(dom_i,"lon"  ,data,lsize) 
    call mct_gGrid_importRAttr(dom_i,"area" ,data,lsize) 
    call mct_gGrid_importRAttr(dom_i,"aream",data,lsize) 
    data(:) = 0.0_R8     
    call mct_gGrid_importRAttr(dom_i,"mask",data,lsize) 
    call mct_gGrid_importRAttr(dom_i,"frac",data,lsize) 
    !
    ! Fill in correct values for domain components
    !
    allocate(work_dom(lsize)) 
    work_dom(:) = 0.0_dbl_kind

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
       do i = ilo, ihi
          n = n+1
          data(n) = TLON(i,j,iblk)*rad_to_deg 
       enddo    !i
       enddo    !j
    enddo       !iblk
    call mct_gGrid_importRattr(dom_i,"lon",data,lsize) 

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
       do i = ilo, ihi
          n = n+1
          data(n) = TLAT(i,j,iblk)*rad_to_deg 
       enddo   !i
       enddo   !j
    enddo      !iblk
    call mct_gGrid_importRattr(dom_i,"lat",data,lsize) 

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
       do i = ilo, ihi
          n = n+1
          data(n) = tarea(i,j,iblk)/(radius*radius)
       enddo   !i
       enddo   !j
    enddo      !iblk
    call mct_gGrid_importRattr(dom_i,"area",data,lsize) 

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
       do i = ilo, ihi
          n = n+1
          data(n) = real(nint(hm(i,j,iblk)),kind=dbl_kind)
       enddo   !i
       enddo   !j
    enddo      !iblk
    call mct_gGrid_importRattr(dom_i,"mask",data,lsize) 

    n=0
    do iblk = 1, nblocks
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
       ihi = this_block%ihi
       jlo = this_block%jlo
       jhi = this_block%jhi
       
       do j = jlo, jhi
       do i = ilo, ihi
          n = n+1
	  if (trim(grid_type) == 'latlon') then
             data(n) = ocn_gridcell_frac(i,j,iblk)
          else
             data(n) = real(nint(hm(i,j,iblk)),kind=dbl_kind)
          end if
       enddo   !i
       enddo   !j
    enddo      !iblk
    call mct_gGrid_importRattr(dom_i,"frac",data,lsize) 

    deallocate(data)
    deallocate(idata)
    deallocate(work_dom)

  end subroutine ice_domain_mct

!=======================================================================
! BOP
!
! !ROUTINE: restart_filename
!
! !INTERFACE:
  character(len=char_len_long) function restart_filename( yr_spec, mon_spec, day_spec, sec_spec )
!
! !DESCRIPTION: Create a filename from a filename specifyer. Interpret filename specifyer 
! string with: 
! %c for case, 
! %t for optional number argument sent into function
! %y for year
! %m for month
! %d for day
! %s for second
! %% for the "%" character
! If the filename specifyer has spaces " ", they will be trimmed out
! of the resulting filename.
!
! !USES:
    use ice_restart, only: runid
!
! !INPUT/OUTPUT PARAMETERS:
  integer         , intent(in)  :: yr_spec         ! Simulation year
  integer         , intent(in)  :: mon_spec        ! Simulation month
  integer         , intent(in)  :: day_spec        ! Simulation day
  integer         , intent(in)  :: sec_spec        ! Seconds into current simulation day
!
! EOP
!
  integer             :: i, n      ! Loop variables
  integer             :: year      ! Simulation year
  integer             :: month     ! Simulation month
  integer             :: day       ! Simulation day
  integer             :: ncsec     ! Seconds into current simulation day
  character(len=char_len_long) :: string    ! Temporary character string 
  character(len=char_len_long) :: format    ! Format character string 
  character(len=char_len_long) :: filename_spec = '%c.cice.r.%y-%m-%d-%s' ! ice restarts

  !-----------------------------------------------------------------
  ! Determine year, month, day and sec to put in filename
  !-----------------------------------------------------------------

  if ( len_trim(filename_spec) == 0 )then
     call shr_sys_abort ('restart_filename: filename specifier is empty')
  end if
  if ( index(trim(filename_spec)," ") /= 0 )then
     call shr_sys_abort ('restart_filename: filename specifier can not contain a space:'//trim(filename_spec))
  end if

  year  = yr_spec
  month = mon_spec
  day   = day_spec
  ncsec = sec_spec

  ! Go through each character in the filename specifyer and interpret if special string

  i = 1
  restart_filename = ''
  do while ( i <= len_trim(filename_spec) )
     if ( filename_spec(i:i) == "%" )then
        i = i + 1
        select case( filename_spec(i:i) )
        case( 'c' )   ! runid
           string = trim(runid)
        case( 'y' )   ! year
           if ( year > 99999   ) then
              format = '(i6.6)'
           else if ( year > 9999    ) then
              format = '(i5.5)'
           else
              format = '(i4.4)'
           end if
           write(string,format) year
        case( 'm' )   ! month
           write(string,'(i2.2)') month
        case( 'd' )   ! day
           write(string,'(i2.2)') day
        case( 's' )   ! second
           write(string,'(i5.5)') ncsec
        case( '%' )   ! percent character
           string = "%"
        case default
           call shr_sys_abort ('restart_filename: Invalid expansion character: '//filename_spec(i:i))
        end select
     else
        n = index( filename_spec(i:), "%" )
        if ( n == 0 ) n = len_trim( filename_spec(i:) ) + 1
        if ( n == 0 ) exit 
        string = filename_spec(i:n+i-2)
        i = n + i - 2
     end if
     if ( len_trim(restart_filename) == 0 )then
        restart_filename = trim(string)
     else
        if ( (len_trim(restart_filename)+len_trim(string)) >= char_len_long )then
           call shr_sys_abort ('restart_filename Resultant filename too long')
        end if
        restart_filename = trim(restart_filename) // trim(string)
     end if
     i = i + 1
  end do
  if ( len_trim(restart_filename) == 0 )then
     call shr_sys_abort ('restart_filename: Resulting filename is empty')
  end if

end function restart_filename

end module ice_comp_mct

